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ABSTRACT 

We present a new formula to numerically construct configurations in rotational equilibrium, which con- 
sist of multiple layers. Each layer rotates uniformly or differentially according to cylindrical rotation- 
laws that are different from layer to layer. Assuming a different barotropic equation of state (EOS) 
for each layer, we solve the Bernoulli equation in each layer separately and combine the solutions by 
^. ■ imposing continuity of the pressure at each boundary of the layers. It is confirmed that a single con- 

tinuous barotropic EOS is incompatible with the junction condition. Identifying appropriate variables 
to be solved, we construct a convergent iteration scheme. For demonstration, we obtain two-layered 
configurations, each layer of which rotates rapidly with either an "f2-constant law" or a "j-constant 
law" or a " u-constant law" . Other rotation laws and/or a larger number of layers can be treated sim- 
ilarly. We hope that this formula will be useful in studying the stellar evolution in multi-dimension 
with the non-spherical configuration induced by rotation being fully taken into account. 
f±j , Subject headings: stars: rotation, stars: evolution, stars: massive 

CO 

^ . 1. INTRODUCTION 

It is well known that stars are generically rotating on the main sequence and, in particular, massive stars are rapid 
O ■ rotators (JFukud a 1982; Tassoul 2000). Although the distribution of angular momentum in the stellar interior is poorly 



h 

star evolves and the central part of the star contracts. In fact some recent theoretical studies on the evolution of 



known except for the sun, it is expected that the inner portion is rotating more rapidly than the outer part as the 



rotating stars have demonstrated that massive s tars in their late evolutionary phases develop a central cor e that is 
rotating more rapidly than the outer envelopes (jHeger et al.ll2000at iHirschi et al.l 120041 iLimongi et al.|[2000l) . Hence 

— , the differential rotation is supposed to exist quite commonly in the stellar interior especially at the advanced stages of 
^- ' evolution. 

\& \ The above-mentioned works on the evolution of rotating stars ignore non-spherical deformations of stars and the 

CO . angle-averaged centrifugal force is added as a correction to the spherical models. This is not justified, however, if the 
star is rotating rapidly and the physical conditions on the rotation axis and on the equatorial plane are substantially 
different. Then the rotational equilibrium should be properly taken into account. This will be particularly important 

\fs for the investigation of the progenitors of gam ma ray bursts, since t hey are supposed to be driven by the gravitational 

(^ • collapse of very rapidly rotating massive stars (jWoosley et al.l 120061 ). 

Over the years substantial effort has been made to numeri cally obtain confi g uratio ns in rotational equilibrium 
in various context s. Be ginning with the pioneering works by lOstriker fe Mark! (|1968l ). a robust iterating for mula 
was developed by [Hachisu (198 6J) and was exten ded to general re lativi stic and/or magnetized stars (iBocquet et al.l 
119951: ICook et all[l994al iKiuchi fc Yoshidall2008bt iKomatsu et"ai][l 989: To mimura fc Erig uchi 20 QJ). The polytropic 



^ equation of state (EOS) was replaced by more generic ones (jCook et al.lll994bi : IKiuchi fc Kotakeil2008at IKiuchi et al.l 
|2009| ). One of the limitations of these studies that hamper the application to the study on the evolution of rotating 
stars, putting aside the treatment of convections and meridional circulations, is the assumption that the whole star 
is rotating cylindrically, that is, the angular velocity is constant on each of th e concentric cy linders and, as a result, 
the EOS is barotropic, that is, the pressure is a function of density alone (|Tassoull |2000[) . On the contrary, the 
theoretical studies based on the spherical models indicate clearly that the stellar core and envelopes composed of 
different elements rotate rather independently of each other, since the steep gradient of mean molecular weight tends 
to suppress the transport of angular momentum beyond the boundary of elements. If one attempts to employ the 
rotational equilibrium configurations in the study of the post-main-sequence evolutions of massive stars, therefore, it 
is almost mandatory to treat multiple layers that obey different rot at ion- laws. 

Motivated by these facts, we present a new formulation to numerically construct non-relativistic configurations in 
rotational equilibrium, which consist of multiple layers. We assume that each layer rotates still cylindrically but the 
rotation-law, namely the angular velocity as a function of the distance from the rotation axis, can be different from 
layer to layer. This assumption allows us to employ the conventional formula based on the Bernoulli equation in 
each layer, a big advantage over the original partial differential equations. The EOS should be barotropic in each 
layer accordingly. We then introduce a junction condition, that is, the continuity of pressure at the boundary of the 
layers to combine them in such a way that the whole star is in rotational equilibrium. In this paper, mainly for 
demonstration purposes, we obtain two-layered configurations with each layer having different polytropic EOS's for 
three representative rotation-laws:(l) J7-constant law (rigid rotation), (2) j-constant law and (3) ^-constant law (see 
the next section for the exact definitions of these rotation laws) . In principle, there is no problem in treating different 
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rotation-laws and/or a larger number of layers. 

The paper is organized as follows. In Sec. [21 we describe the new formulation to obtain rotational equilibrium 
configurations with multiple layers, explaining numerical issues in detail. Section [3] is devoted to the demonstration of 
some model computations. In Sec. |4l we give some discussions and summarize the paper. 

2. BASIC EQUATIONS & ITERATION SCHEME 

2.1. Formulations 

In this paper the configurations in rotational equilibrium are assumed to be axisymmetric and steady with dissipative 
processes being neglected. Possible meridian flows are also ignored. Then the basic equations are 

—^— V l p(r,9) = -\7^ g (r,9)+rsin9n 2 (r 1 9)V l {rsin9), (1) 

V 2 cj> g (r,e) = 4nGp(r,9), (2) 

where the spherical coordinates are used and the subscript i refers to each component. The gravitational potential, 
mass density, pressure, angular velocity and gravitational constant are denoted by <p g , p, p, fl and 67, respectively. For 
a barotropic EOS, i.e., p — p(p), which we assume for each layer throughout this paper, the left hand side of Eq. (JTJ) 
can be integrated. The second term on the right hand side (RHS) of Eq. (TTJ, on the other hand, can be also integrated 
if the angular velocity is a function of the cylindrical radius, i.e., 51 = $7(rsin#). As representative cases, the following 
three rotation laws are chosen in this paper although there is no limitation in principle: 

Q(rsm8) = Hq (rigid rotation), (3) 

Q(r sin 9) — — ^ (constant specific angular momentum in the limit of A — ¥ 0), (4) 

r 2 sin 9 + A 2 

vg 

(r 2 sin 2 9 + A 2 ) 1 / 2 

In these expressions, O , jo and vq are constants that specify how fast the rotation is for each rotation law whereas 
A is a constant that gives a radius of the cylinder, inside which the rotation is almost rigid. The second and third 
rotation laws are referred to as the j-constant and w-constant laws, respectively, according to their limits of A — > 0. 
The first rotation law is also called il-constant law. As already mentioned, the constants Oo, jo, vq and A can be 
different from layer to layer for multi-layered configurations. The values of these constants in the i-th layer (see Fig.[T]) 
are represented by the subscript (i) in the following. 
The integration of Eq. ([1]) in each layer gives the so-called Bernoulli equation, 

H (i)(p( r ,°)) = -<M r >0) + h l(i)^ot{r sin 9 : A {i) ) + c {i} , (6) 

where Hu\ is a specific enthalpy defined by J dP '/ ' p as a function of the mass density alone and cu\ is called a Bernoulli 
constant. The second term on RHS corresponds to the rotational potential and is given for each rotation law as 

h (i) = ^o(i). 9W(rsin0) = -r 2 sin 2 9 (rigid rotation), (7) 



fi(r sin 9) = r ^ — (constant rotational velocity in the limit of A — ► 0). (5) 



hoti) =3o(i), 4>Tot(rsin6 : A {l) ) = - : — — (j-constant law), (8) 

2(r 2 sm 9 + A z ' 

1 



(i)> 



ho(i) = «o(t) > <£rot (r sin 9 : A (l) ) = ^ ln(r 2 sin 2 9 + A 2 t) ) (v-constant law) , (9) 

where the amplitude, Hquj, is expressed separately for later convenience. Note that the Bernoulli constants can be 
different from layer to layer. It is well known that the barotropic condition is equi valent to the r equirement that the 
angular velocity be a function of the cylindrical radius alone, i.e., £1 = Q(rsin9) (|Tassoulll2000f) . In this article, we 
assume a different rotation law for each layer, which means that the angular velocity is not a function of the cylindrical 
radius alone even if the rotation law is cylindrical in each layer. As a result, the EOS cannot be a single continuous 
barotropic one. In fact, the angular velocity is discontinuous across the layer boundary, which leads in general to the 
discontinuity in density as shown later. On the other hand, the pressure is continuous at the layer boundary, which 
can be understood as follows: Multiplied by the density, Eq. ((T|) is written at the layer boundary as 

S7 lP = -pVi<f)g + pr sin 9 ^ 2 4) V,(r sin 9). (10) 

The right hand side (RHS) of this equation contains step functions, that is, the density and angular velocity. Note 
that the gravitational potential, which is obtained by the integration of the density, is continuous. Then the pressure 
is also continuous because otherwise the left hand side of Eq. (fT0|) would give a delta function. This in turn leads to 
the conclusion that a single continuous barotropic EOS cannot be applied across the layer boundary, since the pressure 
could not be continuous for the discontinuous density for such an EOS. Therefore, EOS's that are barotropic in each 
layer but different from layer to layer are required. As the simplest example, we employ in this article polytropic EOS's 
with a different polytropic constant and/or index in each layer. The essential points of our formula are summarized 
as follows: (1) the rotational equilibrium is locally ensured by the Bernoulli equation in each layer and (2) the layer 
boundary is the location where the layers are joined so that the pressure should become continuous. 
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2.2. HSCF scheme for single-layered configurations 

The problem is now reduced to the solution of Eqs. © and © and the search of the location where the pressures 
of the different layers coincide. Before discussin g multi-layered configurations, we briefly review the Hachisu Self- 
Consistent Field (HSCF) scheme (jHachisu l ll986f ), which is known to be a very robust algorithm to solve iteratively 
Eqs. © and ([5} for single-layered configurations in rotational equilibrium and on which our formula is based. In this 
scheme, we first introduce the following non-dimensional variables: 

P = p/pmax, (11) 

f = r/r e , (12) 

ft-o^ot = /io^ot/47rGp max rg, (13) 

c = c/47rGp max r2, (14) 

<j) g = g /47rGp max rg, (15) 

P = p/Pma,xC 2 , (16) 

H = H/c 2 , (17) 

where p max , Pmux and c are the maximum density, pressure and speed of light, respectively, and the subscript i is 
dropped in ho, A, c and H. The radius is normalized by the equatorial radius of the equilibrium configuration, r e , 

which is unknown a priori and is expressed as r e — \ ^^ — by the introduction of a new variable fj. Then 

V P 47rGp^ ax 

Eqs. ([6]) and ([2]) are reduced to 



H{p{f, 6)) = -<t> g {r, 9) + h z fa ot (r sinfl : A) + c, (18) 

Max 

v 2 4> g {f,e) = p{?,e). (19) 

In Eqs. (J18I) and (|19[) . we have two unknown functions <f> g (f,9) and p(f,9) and three constants fj, ho, and c once the 
EOS and rotation law ([7][9]) are specified. It is noted that the specific enthalpy H is a function of the density alone 
because of the barotropic condition. In the HSCF method, we give p max , the equatorial radius f e and polar radius f p 

instead of fj, ho, and c to specify the model and the latter three are treated as unknown variables to be solved. Note 
that f e is unity by the definition of fj (see Eq. (fT2")l ). This choice of variables is essential for the HSCF scheme. Indeed 

other choices such as p max and ho (and f e — 1) fail to obtain convergence in the iteration (see below) more often than 
not. 

The two unknown functions (f> g (f,9) and p(r,9) and three unknown constants fj, ho, and c are obtained iteratively 
in the HSCF method as follows. First we give a trial density distribution p and solve Eq. (fj"9j) to obtain <f> g . As a 
second step, Eq. (|18|) is evaluated at the following points: the stellar surfaces on the equator and on the rotation axis 
and the point of the maximum density, which are denoted as E, P and C ' , respectively (see Fig. [1]). 

(E) = -4> g (l,Tr/2) + h 2 4> rot (l:A) + c, (20) 

(P) Q=-$ g (r p ,0)+h 2 J Iot (0:A)+c, (21) 

(C) -J-H(p max ) = -4> g (fc,n/2) + h 2 4> Iot (fc:A)+c, (22) 

Pmax 

where we made use of the fact that the enthalpy vanishes on the stellar surface. Note that the radius, fc, of point C 
is not known a priori. For the rigid rotation ro t is independent of A. We solve Eqs. (|20|) and ([2~T|) with respect to c 

and ho for 4> g (r, 0) obtained in the first step. With these values of ho and c, we then search for the location where the 
RHS of Eq. (1221 takes the maximum value. The maximum of the RHS of Eq. (|2"2"]l thus obtained gives fj in turn. We 

are now in a position to update p(f,6), solving the Bernoulli equation (|18[) for (f> g (f,9),fJ,ho, and c obtained so far. 
We then repeat the procedure until a sufficient convergence is achieved. 

2.3. Extension to multi-layered configurations 

2.3.1. choice of variables 

We move on to the multi-layered case. Although, for simplicity, we consider only two-layered structures in the 
following, the extension to configurations with a larger number of layers is straightforward. Then we have again two 
unknown functions cf> g (f,9) and p(f,8) and this time five constants fj, h (i), and c^\ in Eqs. © and ©. Another 
important function to be determined is the layer boundary expressed by a function fb(9). Just as in the single layer 
case, the choice of variables and the iteration scheme are critically important to make the scheme convergent. We 
first write down the equations employed to obtain the unknown constants, which correspond to Eqs. (|20l) - ([22|) for the 
single-layered case. As explained above, the pressure should be continuous across the layer boundary. Employing this 
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condition on the equator (point E\ 2 in Figure [TJ and on the rotation axis (point P\ 2 in the same figure) , we write 
down the Bernoulli equation for both sides of the layer boundary at these points: 

(layer 1 side at E 12 ) H (1) (p(r El2 , tt/2)) = -</> 5 (r Bl2 , tt/2) + hl^<p TO t(rE 12 ■ ^4(i)) +c (1) , (23) 

(layer 2 side at E i2 ) iJ (2 )(p(r_E 12 ,7r/2)) = -<p g (rE 12 , tt/2) + /io(2)0rot(r\E 12 : A (2) ) +C( 2 ), (24) 



/'] 



max 



(layer 1 side at P i2 ) H {1) (p(r Pl2 ,0)) = -<&,(rp 12 ,0) + /io (1) ro t(O : A (1) ) +c (1) , (25) 



P, 



max 



(layer 2 side at P 12 ) H {2) (p(r Pl2 ,0)) = -0 9 (rp 12 ,O) + /lo( 2 )0rot(O : A (2) ) +c (2) , (26) 

where we omit * over the normalized variables for notational simplicity and te 12 , rp 12 are the radii of points E\ 2 and 

Pi 2 , respectively. Note that Hn\(p) ^ H^ 2 )(p) because the EOS's are different from layer to layer. We employ the 
previous conditions at points E, P and C, which are written as 

(E) 0=-^(l ! 7r/2) + / l 2 (1) ^ rQt (l:A (1) )+c (1) , (27) 

(P) O = -0 ff (r p ,O) + ^ (1) rot (O:%)) + c (1) , (28) 

(C) -^—H {l) ( Pmax ) = -cb g (rc,Tr/2) + h 2 0il) 4, ot (rc:A {l) ) + c il) . (i = 1 or 2) (29) 

Pmax 



Note that in Eq. (|29|) we do not know a priori in which layer the maximum density point C lies. Using Eqs. (|23|) - (|29[) 
we can determine for the given EOS's and rotation laws seven unknown constants out of /3, /»o(i)i ^0(2), c (i), c (2), 
Pe 12 , Pp 12 , r p, r E 12 , r Pi 2 i m which pe 12 = p(r-£ 12 ,7r /2) and pp 12 = p(rp 12 ,0). This implies that one can give three 
constants to specify the model. As argued shortly, however, pp 12 and rp 12 as well as pe 12 and rg 12 can not be specified 
independently. This can be understood by considering a non-rotating but twodayered configuration. In this case one 
can construct an equilibrium configuration by integrating Eq. (JXJ) radially from the center to rp 12 with the use of the 
gravitational potential L 4irr' pdr' /r 2 . Then it is obvious that pp 12 depends on rp 12 and vice versa (note that the 
maximum density is also fixed). Going back to the multidayered rotational case, we have found that the combination 
of r p , rp 12 and te 12 is a good choice. The triplet of r p , pp 12 , and pe 12 can be an alternative. Note that the inclusion 
of r p seems to be mandatory as has been demonstrated by the HSCF scheme. To summarize, Eqs. (l23l) - ([29|) are used 
to obtain either j3, C(i), c (2 ), /i (i), ^0(2), Pe 12 and p Pl2 for given r p , r El2 and r Pl2 or /?, c (1) , c (2 ), /io(i), ^0(2), r\E 12 and 
rp 12 for given r p , ps 12 and pp 12 . 

2.4. iteration scheme 

Now we proceed to the iteration scheme proposed in this paper. After solving the Poisson equation for the trial 
density distribution, Eqs. (|23l) - (j2"9")l are solved for the variables chosen in the previous section. The procedure is divided 
into the following three steps: 
(1) h (i) and C(i) are calculated from Eqs. (J2"T)) and (|2"5f as 

,2 g (l,7r/2)-0 g (r p ,O) ,„„, 

-</> g (l,7r/2)0 rot (O : A {1) ) + </> g (?>,0)</>rot(l : A {1) ) 
C(1) ro t(l:A (1) )-^ rot (O:A (1) ) ' ( ' 



(2) For a trial value of (3, pe 12 and pp 12 are obtained from Eqs. (|2"3"]l and (j2"5)l and, combined with Eqs. (f2"4")l and (|2~5j) . 
give /i (2) as 

,2 _ 0rot(rg 12 : A (1) ) - rot (O : A (1) ) u2 
° (2) "^ot(r Bl2 :A (2) )-0 rot (O:A (2) )Ni) 

#(2) (Pgi 2 ) - g(l) (P£?i 2 ) - g(2) (PP 12 ) + ^(1) (PP 12 ) g , 32 > 

</>rot(VE 12 : A( 2 )) — 0rot(O : ^4(2)) Pmax 

Then c (2) is obtained from Eq. flU} or (j2"6l . 

(3) The point that gives the largest value to — g (r, tt/2) + h^^tfirotir : ^4(i)) +C(i) is searched in layer 1 and is referred 
to as point CI. The counter part in layer 2 is then looked for and is called point C2. The point C is found by comparing 
CI and C2 to E\ 2 . The maximum value thus obtained is divided by ff(i)(/9 m ax)/pmax with the appropriate i and the 
updated value of /3 is obtained. The steps (2) and (3) are repeated until the value of (3 converges at a sufficient level. 

2.4.1. Layer boundary 
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Table 1 
Comparison of the non-rotational two-layered configurations obtained by the ordinary method (ID) and 

THE PRESENT FORMULA (2D). THE RADIUS, MASS IN CGS UNIT AND NORMALIZED VIRIAL ARE DENOTED BY R, M, AND VC, 
RESPECTIVELY. FOR THIS MODEL WE ADOPT Te 12 = rp 12 = 1.718 X 1CT 2 . 



R [cm] M [M ] c (1) C( 2 ) VC 

ID 3.218E+08 4.292E+00 -3.622E-03 -8.579E-03 

2D 3.224E+08 4.291E+00 -3.601E-03 -8.533E-03 2.815E-05 



The final step of the iteration is the updates of the density distribution and the location of the layer boundary. This 
is accomplished as follows. Regarding the specific enthalpy as a function of pressure alone, we first solve the Bernoulli 
equation, Eq. (|6|), for each layer in the absence of the other layer and the layer boundary as a result: 

P{1) (r,9)=H- 1 1 ) ((-4,^,9) + hl {1) <p mt (r sine : A {1) ) + c w )^\ , (33) 

p {2) (r,6)=H^ ^ g (r,e) + hl {2) cf>Ur sin9 ■.A (2) ) + c {2) )^^j , (34) 

where H7.-. is an inverse function of the specific enthalpy for layer i. Since the pressure should be continuous at the 
layer boundary as discussed in Sec. I2.1[ we look for a point on each radial ray (a line with 9 = const.), at which pn\ 
and p( 2 ) coincides with each other. This gives the updated layer boundary as r = Tb{0). On the other hand, pn)(r, e) 
and P(2)( r : e) obtained above give the updated density distribution for each layer. This closes the iteration procedure. 
We go back to the Poisson equation and repeat all the steps until a sufficient level of convergence is reached. 

3. RESULT 

To demonstrate that the new formula described above really works, we will apply it to two-layered configurations for 
some representative rotation laws. For simplicity, we employ two polytropic EOS's, Pi(p) — Kip 1+1 ' ni , in which Ki and 
Hi are the polytropic constants and indices for layer i. We take rather arbitrarily K\ — 1.015x 10 15 and K 2 = 5.14x 10 14 
in cgs unit and m = n 2 = 3. Note that the specific enthalpy is given as H^ = (m + l)Kip l / ni = n^if™*'™ 1 pi/i+«. 
for layer i. Adopting p max = 5.629 x 10 9 g/cm 3 we obtain white dwarf-like configurations. 

We work with the normalized variables (Eqs. (|ll[) - (|17p ) and the numerical domain covers a quadrant of the meridian 
section, ^ r ^ 1 and 0^0^ it/2, under the assumption of equatorial symmetry. We typically deploy 1000 mesh 
points on the r-coordinate and 200 grid points on the ^-coordinate to obtain an acceptable accuracy, which is confirmed 
by the normalized virial equation (|Cowl ing 1965) defined as 

VC=\2T + W + 3U\/\W\, (35) 

where T, W, and U are the rotational, gravitational and internal energies, respectively. The ratio should v anish for 
exact solutions. Note finally that the Poisson equation is solved by the Green function method employed bv iHachisul 
(fl986l) . 

3.1. Non-rotational case 

As a mandatory step, we first construct a non-rotational configuration with two layers according to the present 
formula and compare it with the solution obtained by the ordinary and much simpler method, that is, the radial 
integration. Setting r p — 1 and te 12 — r P 12 , we obtain /io(i) = from Eq. ([30]) and pe 12 = Pp± 2 fr° m Eqs. (|23|) and 
(|2"5|) that are actually identical with each other. Then, Eq. (|3"2")l gives /io(2) = because the first term vanishes owing 
to fto(i) = and the second term is also zero because of pe 12 — PP\ 2 ■ As a consequence, Eq. (|24p becomes identical 
with Eq. (121)1) . It should be noted that c^ ^ C( 2 ) even in the non-rotational case, since the EOS's are different between 
two layers. 

Table [l] summarizes the comparison of the non-rotating configuration obtained by the new formulation, which is 
referred to as " 2D" in the table, with the solution of the one-dimensional hydrostatic equation obtained by the radial 
integration from center to surface, which is called " ID" . The radius, mass and the Bernoulli constants in 2D, which 
characterize the equilibrium configuration, nicely agree with the ID counter parts. Indeed, the relative errors in these 
quantities are less than one percent and the normalized virial is of the order of 10~ 5 , both of which indicate that the 
present method can reproduce the non-rotational configuration. 

3.2. Q-constant case 

We now proceed to the rotational cases. In this section we deal with the simplest one, that is, the combination of two 
rigid rotations (£7-constant laws in Eq. (j3|)). We have constructed two configurations with (r p , te 12 , rp 12 ) = (0.9, 0.5, 0.5) 
and (0.8,0.5,0.433). In Table [5J we show some key quantities that characterize the configurations. As expected, the 
angular velocities are different between the layers. Note that in the present formula, the angular velocity is not specified 
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Table 2 
Radius, mass, angular velocities, Bernoulli constants and virial for two £7-constant configurations. 

Qp,r El2 ,rp 12 ) R [cm] M [M©] U m [rad/s] (2 ) [rad/s] c (1) c^g) VC 

(0.9,0.5,0.5) 2.802E+08 1.554E+00 1.441E+00 1.087E+00 -2.207E-03 -2.911E-03 2.120E-05 
(0.8,0.5,0.433) 3.116E+08 1.585E+00 1.854E+00 2.910E+00 -1.838E-03 -2.453E-03 1.499E-05 

Table 3 
Radius, mass, specific angular momenta/rotational velocities, Bernoulli constants, and virial for 

j-constant and u-constant cases. 

j-constant case 

(r P ,r El2 ,rp 12 ,A (1) ,A(2)) R [cm] M [M Q ] j 0(1) [cm 2 /s] j 0(2) [cm 2 /s] c (1) eg VC 

(0.6,0.5,0.32,0.1,0.05) 2.405E+08 1.739E+00 2.683E+16 1.435E+16 -3.501E-03 -4.956E-03 4.969E-05 

(0.6,0.5,0.5,0.1,0.1) 2.555E+08 1.711E+00 2.748E+16 2.074E+16 -2.870E-03 -4.042E-03 4.404E-05 

v-constant case 

(r P ,r El2 ,rp 12 ,A (1) ,A(2)) R [cm] M [M Q ] » 0(1 ) [cm/s] »o( 2 ) [cm/s] c (1) eg VC 

(0.6,0.5,0.32,0.1,0.05) 3.037E+08 1.823E+00 4.694E+08 4.709E+08 -1.837E-03 -2.434E-03 2.204E-05 

(0.6,0.5,0.5,0.1,0.1) 3.057E+08 1.654E+00 4.506E+08 3.412E+08 -1.632E-03 -2.295E-03 1.326E-05 

but solved. It is te 12 and rp 12 that dictate the angular velocities. If one wants to construct a configuration for particular 
angular velocities, another iteration is needed for shooting. Incidentally, we can impose in principle the combinations 
that sat isfy re,, < rp ,^ , which is intuitively unlikely because rotations tend to flatten equilibrium configurations in 
general (|Hachisullf9 86). In fact, we find that the solutions for such cases have a negative centrifugal force, which is, of 
course, unphysical. We thus confirm that the inner layer is still oblate in the multi-layered configurations. 

Figure [2] displays the contour plots of the density and angular velocity in the meridian section for the model with 
{ r pi r E 121 rp 12 ) — (0.8,0.5,0.433). The thick curve in the figure represents the layer boundary. It should be noted 
that the density is discontinuous at the layer boundary as mentioned earlier. This is generally the case. On the 
other hand, the pressure should be continuous across the boundary. To check this, we show in Fig. [3] the density 
and pressure profiles along the radial lines with 9 = 0, 7r/4 and 7r/2. It is clear that the density is discontinuous at 
1.2 x I0 8 cm < r < 1.6 x I0 8 cm whereas the pressure profiles are continuous in the same region. From this we can 
understand again why different EOS's are needed for each layer. 

The values of VC in Tableware of the same orders of magnitude, ~ 1 0~ 5 , as in the non-rotational case. This implies 
that the configurations we have constructed by the new formula are in rotational equilibrium to the same accuracy 
as the spherical configuration in the previous section is in hydrostatic equilibrium. Figure U plots the values of VC 
as a function of the number of grid points. It is clearly demonstrated that the accuracy is increased as the resolution 
becomes better although the convergence is rather slow. We infer that this slow convergence is due to the discontinuity 
in the density distribution, which enters the virial equation through the gravitational binding energy W and rotational 
energy T. The important point, however, is the fact that the accuracy is improved as the grid number increases. We 
are thus confident that our new formulation has indeed succeeded in finding two-layered configurations in rotational 
equilibrium. 

3.3. j-constant and v-constant cases 

We move on to the combinations of other rotation laws. In the following configurations, each layer is rotating 
differentially. Although various combinations are actually possible, only those in the same rotation law are considered 
here just for simplicity. Thus the j-constant case refers to the configurations, in which each layer obeys the j-constant 
law in Eq. (j4|) although ho and A are different between the layers. The same is true of the v-constant case. For 
each case we have constructed two configurations that have (r p ,rE 121 rp 12 , A^, A^ 2 )) — (0.6,0.5,0.32,0.1,0.05) and 
(0.6,0.5,0.5,0.1,0.1). Note again that h ^, or how fast each layer is rotating, are not specified but solved and te 12 
and rp 12 are the control parameters. On the other hand, the degree of differential rotations or Au\ can be specified 
freely. 

Table [3] summarizes the quantities that characterize the configurations whereas Figs. [5] and [6] display the contour plots 
of the density and angular velocity for the models in the j-constant and u-constant cases, respectively. As expected 
and shown in these figures, the differentially rotating models are more deformed than the one presented in the previous 
section, in which each layer rotates rigidly. The values of VC in Table [3] are of the same order, ~ 10~ 5 , as in the 
non-rotational and f2-constant cases, that fact indicates our formulation's capability of constructing configurations 
with strongly differential rotations. This may be important in dealing with the progenitors of GRB ( Wo oslev et al.l 
[2006T ). 

4. SUMMARY AND DISCUSSIONS 

Bearing in mind the application to the study of rotational massive stars in their late evolutionary phase, in this 
paper we have proposed a new formula to construct multi-layered configurations in rotational equilibrium. This is 
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an extension of the Hachisu self-consistent field scheme that is based on the Bernoulli equation and meant originally 
for single-layered configurations that are rotating cylindrically with a barotropic EOS. In our method, on the other 
hand, each layer is assumed to rotate still cylindrically with a barotropic EOS but the rotation laws and EOS's are 
different from layer to layer. We have shown that the pressure should be continuous at the layer boundary whereas the 
density is in general discontinuous across the boundary, which is an alternative demonstration that the EOS cannot 
be identical for the adjacent layers. We have identified the variables that are appropriate to make the iteration scheme 
convergent. This is indeed a crucial ingredient in our formula. 

For demonstration, we have actually constructed several configurations with two layers for three representative 
rotation laws, which we have referred to as the J7-constant, j-constant and u-constant laws in this paper. We have 
found that the virial equation is satisfied with a typical error of 10 -5 irrespective of the rotation laws if we deploy 
1000 x 200 mesh points and we have also demonstrated that the error is reduced as the number of mesh points is 
increased. Incidentally, it has been confirmed that a non-rotational configuration is also reproduced by the present 
scheme. From these results it is obvious that our method works well and is robust indeed. The application of the 
present formula to more realistic problems will be published elsewhere (|Nagakura et al.ll2010[ ). 

As commented in Sec. [2j it is straightforward to extend our scheme to the configurations with more than two layers 
though the procedure becomes a bit more involved. Although we have combined the rotation laws of the same family 
but with different parameters for simplicity in this paper, two rotation laws of different families can be treated in 
the same way. The implementation of more realistic EOSs w i ll pose no problem in p rinciple as long as they are 
barotropic. We may employ the idea by Uackson et al.l (|2005| ): iMacGregor et al.l (|2007| ) that the pressure, density, 
and temperature are assumed to be functions of the effective potential alone. Moreover, the present formula will be 
able to treat configurations with a topology of torus by relaxing the assumption that the surface extends itself to the 
symmetry axis and by choosing appropriate points on the equator to impose the conditions corresponding to Eqs. (1251) . 
(|26|) and ([28| although we do not know how realistic such configurations are. 

In our formula, the layer boundary is determined from the condition that the pressure be continuous there. In 
reality, however, the layers in the stellar interior correspond to the regions of different chemical compositions and their 
boundaries are determined by the thermodynamical conditions for nuclear burnings. This difference originates from 
the fact that we have imposed piece-wise cylindrical rotation laws. In the actual stellar interior, each layer obeys a 
boroclinic EOS and, as a result, rotates non-cylindrically. Moreover, the gas motions in the meridian section such as 
convections and meridional circulations are likely to exist generically. Then the original partial differential equations 
should be solved somehow, which is a formidable task and will need an entirely new approach. Our formula, therefore, 
is admittedly a rather crude approximation to the reality but, hopefully, not so bad one if one chooses an appropriate 
rotation law for each layer. In fact, it will be much better than any approximate configurations with only a single-layer. 

The real challenge will be to somehow implement chemical evolutions to the sequence of rotational configurations. 
One possibility may be an extension of the idea employed in most of the current one dimensional evolution models 
of rotational massive stars (|Heger et al.ll2000at iHirschi et al.l 120041 : iLimongi et al.ll2000l ). Under the assumption that 
the thermodynamical conditions as well as the chemical abundances are uniform on each surface of constant effective 
potential, we solve the nuclear network locally and then transfer generated energy on the multi-dimensional mesh. 
The transport of angular momentum may be also approximated by diffusion. Since the resultant distributions of 
thermodynamical quantities and elements will in general be non-uniform on the surface of constant effective potential, 
we will take their angular averages on the surface and solve the new rotational equilibrium for the obtained equations 
of state and rotation law. This completes the single cycle and the iteration of this process will give the temporal 
evolution of rotational stars. We hope that this procedure is feasible and that the formulation presented in this paper 
will contribute to the study of the influences of non-sphericity on the evolution of rapidly rotating massive stars. 



Numerical computations were in part carried on XT4 and general common use computer system at the center for 
Computational Astrophysics, CfCA, the National Astronomical Observatory of Japan and on NEC-SX8 at Yukawa 
Institute for Theoretical Physics in Kyoto University. This study was supported in part by the Grants-in-Aid for the 
Scientific Research from the Ministry of Education, Science and Culture of Japan (Nos. 80251403 and 19104006). 
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r sin 9 

Fig. 1. — Schematic picture of the two-layered configuration in rotational equilibrium in the meridional section. The solid and dashed 
curves represent the layer boundary and stellar surface, respectively. The circles attached by letters C, E12, E, P12 and P correspond 
to the point of density maximum, layer boundary and surface on equator and rotation axis, respectively. Each layer is referred to by the 
number 1 or 2. See the text for more details. 
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Fig. 2. — Contour plots in the meridional section of the density (left panel) and angular velocity (right panel) for the J7-constant case. 
The contour levels are equally spaced in the logarithmic scale in panel (a) and in the linear scale in (b). The thick curves in the stellar 
interiors correspond to the layer boundary in both panels. For this model (t p ,te 12 >' r Pi2) = (0.8,0.5,0.433) are adopted. Each layer is 
rotating rigidly. 
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Fig. 3. — Density (left panel) and pressure (right panel) profiles for the model in Fig. [2] In both panels, the solid, dashed and dotted 
curves show the distributions along the radial lines with 8 = 0, n/A and 7r/2, respectively. The insets are the magnification of the vicinity 
of layer boundary. 
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Fig. 4. — Convergence test. The normalized virial is shown as a function of the number of grid points for the f2-constant model with 
(7>! r £i2! r -Pi2) = (0.8,0.5,0.433). N r and N$ denote the numbers of grid points on r- and ^-coordinates, respectively. 
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Fig. 5. — Contour plots in the meridional section of the density (left panel) and angular velocity (right panel) for the jf-constant case 
with (r p , te 12 ! r P 12 ) = (0.6, 0.5, 0.32), A(i) = 0.1 and A(2) = 0.05. The contour levels are equally spaced in the logarithmic scale in both 
panels. The thick curves in the stellar interiors correspond to the layer boundary in both panels. 
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Fig. 6. — Contour plots in the meridional section of the density (left panel) and angular velocity (right panel) for the ^-constant case 
with (r p ,rE 12 , rp 12 ) = (0.6, 0.5, 0.5), An\ = 0.1 and An) = 0.1. The contour levels are equally spaced in the logarithmic scale in both 
panels. The thick curves in the stellar interiors correspond to the layer boundary in both panels. 



